######################################################################
# Politics or Performance? Leadership Accountability in UN Peacekeeping
# by Magnus Lundgren, Kseniya Oksamytna, and Vincenzo Bove
# Replication script, R
# June 4, 2021
######################################################################


# Load required packages; use "install.packages" function as necessary
library(survival)
library(stargazer)
library(gginnards)
library(survminer)
library(simPH)
library(lme4)


# Load data 

dataFC<-read.csv("~/Dropbox/Tenure paper/Supplementary online material/dataFC_repl.csv")
dataSRSG<-read.csv("~/Dropbox/Tenure paper/Supplementary online material/dataSRSG_repl.csv")


# Create survival objects

survobjectFC<-Surv(dataFC$FCtenure, dataFC$FCtenure2, dataFC$FCexit)
survobjectSRSG<-Surv(dataSRSG$SRSGtenure, dataSRSG$SRSGtenure2, dataSRSG$SRSGexit)


# FIGURES AND TABLES IN ARTICLE

# Figure 1

temp<-aggregate(FCduration~ForceCommander,data=dataFC[dataFC$FCrightcensored==0,],FUN=mean)
temp$FCrightcensored<-0
temp2<-aggregate(FCduration~ForceCommander,data=dataFC[dataFC$FCrightcensored==1,],FUN=mean)
temp2$FCrightcensored<-1
temp<-rbind(temp,temp2)
temp3<-aggregate(SRSGduration~SRSG,data=dataSRSG[dataSRSG$SRSGrightcensored==0,],FUN=mean)
temp3$SRSGrightcensored<-0
temp4<-aggregate(SRSGduration~SRSG,data=dataSRSG[dataSRSG$SRSGrightcensored==1,],FUN=mean)
temp4$SRSGrightcensored<-1
temp3<-rbind(temp3,temp4)
par(mfrow=c(1,2))
hist(temp$FCduration,col="grey60",main="Force commanders",ylab="Event count",xlab="Tenure in months",xlim=c(0,70),ylim=c(0,35),breaks=20)
hist(temp$FCduration[temp$FCrightcensored==1],add=T,col="grey")
legend("topright",c("Exit","Right-censoring"),fill=c("grey60","grey"),cex=.8,bty="n")
hist(temp3$SRSGduration,col="grey60",main="SRSGs",ylab="Event count",xlab="Tenure in months",xlim=c(0,70),ylim=c(0,35),breaks=20)
hist(temp3$SRSGduration[temp3$SRSGrightcensored==1],add=T,col="grey")
legend("topright",c("Exit","Right-censoring"),fill=c("grey60","grey"),cex=.8,bty="n")


# Figure 2

keeps<-c("SRSGtenure","SRSGtenure2","SRSGexit")
data<-dataSRSG[keeps]
colnames(data)<-c("tenure","tenure2","exit")
data$type<-"SRSG"
keeps<-c("FCtenure","FCtenure2","FCexit")
data2<-dataFC[keeps]
colnames(data2)<-c("tenure","tenure2","exit")
data2$type<-"FC"
data<-rbind(data,data2)
survobject<-Surv(data$tenure, data$tenure2, data$exit)
fit<-survfit(survobject~data$type)
par(mfrow=c(1,1))
plot(fit, xlim=c(0,60), ylab="Probability of leader remaining in post", xlab="Tenure in months", lwd=2, main="",lty=c(1,3)) # how take out confidence intervals
legend("topright",c("FCs","SRSGs"),lty=c(1,3),cex=.8,lwd=c(2,2))


# Table 1

mod1<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T)

mod2<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3:FC_P5+cluster(missionID), data = dataFC,robust=T)

mod3<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3:log(gle_gdp)+cluster(missionID), data = dataFC,robust=T)

mod4<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3:Troops_log+cluster(missionID), data = dataFC,robust=T)

mod5<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataSRSG,robust=T)

mod6<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3:SRSG_P5+cluster(missionID), data = dataSRSG,robust=T)

mod7<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3:log(gle_gdp)+cluster(missionID), data = dataSRSG,robust=T)

mod8<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3:Troops_log+cluster(missionID), data = dataSRSG,robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]
robust3     <- summary(mod3)$coefficients[,4]
robust4     <- summary(mod4)$coefficients[,4]
robust5     <- summary(mod5)$coefficients[,4]
robust6     <- summary(mod6)$coefficients[,4]
robust7     <- summary(mod7)$coefficients[,4]
robust8     <- summary(mod8)$coefficients[,4]

stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,digits=3,se=list(robust1,robust2,robust3,robust4,robust5,robust6,robust7,robust8),type="text",df=F, keep.stat=c("n","ll","lr"))


# Figure 3

mod1<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T)
Sim1 <- coxsimLinear(mod1, b = "performance3", qi="First Difference", Xj = seq(0, 100, by = 1),  nsim = 1000,ci=.95) 
f1<-simGG(Sim1, xlab="Mission performance index score (FCs)",ylab="Percentage change in hazard rate",lcolour="black",pcolour="grey50")
f1+geom_rug(alpha=.5,position="jitter",sides="b")
which_layers(f1,"GeomRibbon")
p<-delete_layers(f1,"GeomRibbon")
p<-p+ geom_line(aes(y = Min), color = "black", linetype = "dashed")
p<-p+ geom_line(aes(y = Max), color = "black", linetype = "dashed") +geom_rug(sides="b",alpha = 1, position = "jitter")
p


# Figure 4

mod1<-coxph(formula = survobjectFC ~  performance3+FC_P5+gle_gdp_log+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T)
new_data<-data.frame("performance3"=mean(dataFC$performance3),"Troops_log"=mean(dataFC$Troops_log),FC_P5=c(0,1), gle_gdp_log=mean(na.omit(dataFC$gle_gdp_log)),FCpreviousexperiencedummy=0,Chapter.VII=0,Troops_log=na.omit(mean(dataFC$Troops_log)),missionend=0,acting=0)
survfit(mod1,newdata=new_data,conf.int=.95)
plot(survfit(mod1,newdata=new_data,conf.int=.90),conf.int=F,lty=c(1,3), xmax=30,lwd=2,main="",xlab="Tenure (months)",ylab="Probability of survival")
legend("topright",lty=c(3,1),legend=c("P5 FCs","Non-P5 FCs"),lwd=2,cex=.7)


# Table 2

mod9<-coxph(formula = survobjectFC ~  CumulativeFatalities_log+CumulativeMalicious_log+cumulative_onesided_best_log+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T)

mod10<-coxph(formula = survobjectFC ~  cum_sea_mil+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+cluster(missionID), data = dataFC,robust=T,control=coxph.control(eps=0.002))

mod11<-coxph(formula = survobjectSRSG ~  CumulativeFatalities_log+CumulativeMalicious_log+cumulative_onesided_best_log+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataSRSG,robust=T)

mod12<-coxph(formula = survobjectSRSG ~  cum_sea_mil+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII++missionend+cluster(missionID), data = dataSRSG,robust=T,control=coxph.control(eps=0.002))

robust9     <- summary(mod9)$coefficients[,4]
robust10     <- summary(mod10)$coefficients[,4]
robust11    <- summary(mod11)$coefficients[,4]
robust12     <- summary(mod12)$coefficients[,4]

stargazer(mod9,mod10,mod11,mod12,digits=3,se=list(robust9,robust10,robust11,robust12),type="text",df=F,keep.stat=c("n","ll","lr"))


# Figure 5

mod1<-coxph(formula = survobjectFC ~  CumulativeFatalities_log+CumulativeMalicious_log+cumulative_onesided_best_log+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T)
Sim1 <- coxsimLinear(mod1, b = "CumulativeFatalities_log", qi="First Difference", Xj = seq(0,2.5, by = .1),  nsim = 1000,ci=.95) 
f2<-simGG(Sim1, xlab="Battle deaths during FC tenure (1,000s) (logged)",ylab="Percentage change in hazard rate",lcolour="black",pcolour="grey50")
f2+geom_rug(alpha=.5,position="jitter",sides="b")
p<-delete_layers(f2,"GeomRibbon")
p<-p+ geom_line(aes(y = Min), color = "black", linetype = "dashed")
p<-p+ geom_line(aes(y = Max), color = "black", linetype = "dashed") +geom_rug(sides="b",alpha = 1, position = "jitter")
p


# FIGURES AND TABLES IN ONLINE APPENDIX

# Table A1

test<-dataFC[c("UNacronym","ForceCommander")]
test<-unique(test)
test$count<-1
x<-aggregate(count~UNacronym,data=test,FUN=sum)
mini<-aggregate(year~UNacronym,data=dataFC,FUN=min)
maxi<-aggregate(year~UNacronym,data=dataFC,FUN=max)
x<-merge(x,mini,by=c("UNacronym"))
x<-merge(x,maxi,by=c("UNacronym"))
test<-dataSRSG[c("UNacronym","SRSG")]
test<-unique(test)
test$count<-1
z<-aggregate(count~UNacronym,data=test,FUN=sum)
x<-merge(x,z,by=c("UNacronym"),all.x=T)
x<-x[,c("UNacronym","count.y","count.x","year.x","year.y")]
colnames(x)<-c("Mission","SRSGs","FCs","Earliest year","Latest year")
x


# Table A2

datasummary<-subset(dataFC, select=c("performance3","cum_good_combined_new_3","cum_bad_combined_new_3","FC_P5","gle_gdp_log","Troops_log","FCpreviousexperiencedummy","Chapter.VII","acting","missionend","DE","vdem_libdem","ctroops_log", "gle_rgdpc_log"))
stargazer(datasummary, digits=2,type="text",covariate.labels = c("Performance (index)","Performance (good, 3 months)","Performance (bad, 3 months)","P5","GDP (log)","Mission troop contribution (log)","Previous experience","Chapter VII","Acting","End of mission", "Diplomatic network","Liberal democracy index","General troop contribution (log)","GDP per capita (log)"))

datasummary<-subset(dataSRSG, select=c("performance3","cum_good_combined_new_3","cum_bad_combined_new_3","SRSG_P5","gle_gdp_log","Troops_log","SRSGpreviousexperiencedummy","Chapter.VII","acting","missionend","DE","vdem_libdem","ctroops_log", "gle_rgdpc_log"))
stargazer(datasummary, digits=2,type="text",covariate.labels = c("Performance (index)","Performance (good, 3 months)","Performance (bad, 3 months)","P5","GDP (log)","Mission troop contribution (log)","Previous experience","Chapter VII","Acting","End of mission", "Diplomatic network","Liberal democracy index","General troop contribution (log)","GDP per capita (log)"))


# Table A3

datasummary<-subset(dataFC, select=c("performance3","cum_good_combined_new_3","cum_bad_combined_new_3","FC_P5","gle_gdp_log","Troops_log","FCpreviousexperiencedummy","Chapter.VII","acting","missionend"))
correlation.matrix <- cor(datasummary,use="pairwise.complete.obs")
stargazer(correlation.matrix, title="Correlation Matrix (FC data)", digits=2,type="text")

datasummary<-subset(dataSRSG, select=c("performance3","cum_good_combined_new_3","cum_bad_combined_new_3","SRSG_P5","gle_gdp_log","Troops_log","SRSGpreviousexperiencedummy","Chapter.VII","acting","missionend"))
correlation.matrix <- cor(datasummary,use="pairwise.complete.obs")
stargazer(correlation.matrix, title="SRSG data", digits=2,type="text")


# Table A4

mod1<-coxph(formula = survobjectFC ~  cum_good_combined_new_3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+cluster(missionID), data = dataFC,robust=T)

mod2<-coxph(formula = survobjectFC ~  cum_bad_combined_new_3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+cluster(missionID), data = dataFC,robust=T,control=coxph.control(iter.max=10000))

mod3<-coxph(formula = survobjectSRSG ~  cum_good_combined_new_3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+cluster(missionID), data = dataSRSG,robust=T,control=coxph.control(eps=0.003,iter.max=10000))

mod4<-coxph(formula = survobjectSRSG ~  cum_bad_combined_new_3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+cluster(missionID), data = dataSRSG,robust=T,control=coxph.control(eps=0.003,iter.max=10000),init=c(-17,.5,0,0,0,-1))

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]
robust3     <- summary(mod3)$coefficients[,4]
robust4     <- summary(mod4)$coefficients[,4]

stargazer(mod1,mod2,mod3,mod4,digits=3,se=list(robust1,robust2,robust3,robust4),type="text",df=F, keep.stat=c("n","ll","lr"))


# Table A5

survobjectFC1989<-Surv(dataFC$FCtenure[dataFC$year>1989], dataFC$FCtenure2[dataFC$year>1989], dataFC$FCexit[dataFC$year>1989])
survobjectSRSG1989<-Surv(dataSRSG$SRSGtenure[dataSRSG$year>1989], dataSRSG$SRSGtenure2[dataSRSG$year>1989], dataSRSG$SRSGexit[dataSRSG$year>1989])

mod1<-coxph(formula = survobjectFC1989 ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC[dataFC$year>1989,],robust=T)

mod2<-coxph(formula = survobjectSRSG1989 ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataSRSG[dataSRSG$year>1989,],robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]

stargazer(mod1,mod2,digits=3,se=list(robust1,robust2),type="text",df=F, keep.stat=c("n","ll","lr"))


# Table A6

mod1<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+ctroops_log+cluster(missionID), data = dataFC,robust=T)

mod2<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+ctroops_log+cluster(missionID), data = dataSRSG,robust=T)

mod3<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+DE+cluster(missionID), data = dataFC,robust=T)

mod4<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+DE+cluster(missionID), data = dataSRSG,robust=T)

mod5<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+log(gle_rgdpc)+cluster(missionID), data = dataFC,robust=T)

mod6<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+log(gle_rgdpc)+cluster(missionID), data = dataSRSG,robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]
robust3     <- summary(mod3)$coefficients[,4]
robust4     <- summary(mod4)$coefficients[,4]
robust5     <- summary(mod5)$coefficients[,4]
robust6     <- summary(mod6)$coefficients[,4]

stargazer(mod1,mod2,mod3,mod4,mod5,mod6,digits=3,se=list(robust1,robust2,robust3,robust4,robust5,robust6),type="text",df=F,keep.stat=c("n","ll","lr"))



# Table A7

survobjectFC_not_last_6<-Surv(dataFC$FCtenure[dataFC$last_6==0], dataFC$FCtenure2[dataFC$last_6==0], dataFC$FCexit[dataFC$last_6==0])
survobjectSRSG_not_last_6<-Surv(dataSRSG$SRSGtenure[dataSRSG$last_6==0], dataSRSG$SRSGtenure2[dataSRSG$last_6==0], dataSRSG$SRSGexit[dataSRSG$last_6==0])

mod1<-coxph(formula = survobjectFC_not_last_6 ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+acting+cluster(missionID), data = dataFC[dataFC$last_6==0,],robust=T)

mod2<-coxph(formula = survobjectSRSG_not_last_6 ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+acting+cluster(missionID), data = dataSRSG[dataSRSG$last_6==0,],robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]

stargazer(mod1,mod2,digits=2,se=list(robust1,robust2),type="text",df=F,keep.stat=c("n","ll","lr"))


# Table A8

tempFC<-dataFC[dataFC$OBSERVER==0,] 
tempSRSG<-dataSRSG[dataSRSG$OBSERVER==0,]

survobjectFC_temp<-Surv(tempFC$FCtenure, tempFC$FCtenure2, tempFC$FCexit)
survobjectSRSG_temp<-Surv(tempSRSG$SRSGtenure, tempSRSG$SRSGtenure2, tempSRSG$SRSGexit)

mod1<-coxph(formula = survobjectFC_temp ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempFC,robust=T)

mod2<-coxph(formula = survobjectSRSG_temp ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempSRSG,robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]

stargazer(mod1,mod2,digits=2,se=list(robust1,robust2),type="text",df=F,keep.stat=c("n","ll","lr"))


# Table A9

tempFC<-dataFC[dataFC$performance3>31.1,]
tempSRSG<-dataSRSG[dataSRSG$performance3>17,]
survobjectFC_temp<-Surv(tempFC$FCtenure, tempFC$FCtenure2, tempFC$FCexit)
survobjectSRSG_temp<-Surv(tempSRSG$SRSGtenure, tempSRSG$SRSGtenure2, tempSRSG$SRSGexit)

mod1<-coxph(formula = survobjectFC_temp ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempFC,robust=T)

mod2<-coxph(formula = survobjectSRSG_temp ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempSRSG,robust=T)

tempFC<-dataFC[dataFC$performance3<max(dataFC$performance3),] 
tempSRSG<-dataSRSG[dataSRSG$performance3<max(dataSRSG$performance3),]

survobjectFC_temp<-Surv(tempFC$FCtenure, tempFC$FCtenure2, tempFC$FCexit)
survobjectSRSG_temp<-Surv(tempSRSG$SRSGtenure, tempSRSG$SRSGtenure2, tempSRSG$SRSGexit)

mod3<-coxph(formula = survobjectFC_temp ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempFC,robust=T)

mod4<-coxph(formula = survobjectSRSG_temp ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempSRSG,robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]
robust3     <- summary(mod3)$coefficients[,4]
robust4     <- summary(mod4)$coefficients[,4]

stargazer(mod1,mod2,mod3,mod4,digits=2,se=list(robust1,robust2,robust3,robust4),type="text",df=F,keep.stat=c("n","ll","lr"))


# Table A10

tempFC<-dataFC[c("FCtenure","FCtenure2","FCexit","performance3","Troops_log","FC_P5","gle_gdp","FCpreviousexperiencedummy","Chapter.VII","missionID","missionend","acting")]
tempFC<-na.omit(tempFC)

tempSRSG<-dataSRSG[c("SRSGtenure","SRSGtenure2","SRSGexit","performance3","Troops_log","SRSG_P5","gle_gdp","SRSGpreviousexperiencedummy","Chapter.VII","missionID","missionend","acting")]
tempSRSG<-na.omit(tempSRSG)

survobjectFC_temp<-Surv(tempFC$FCtenure, tempFC$FCtenure2, tempFC$FCexit)
survobjectSRSG_temp<-Surv(tempSRSG$SRSGtenure, tempSRSG$SRSGtenure2, tempSRSG$SRSGexit)

mod1<-coxph(formula = survobjectFC_temp ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempFC,robust=T)
dfbeta <- residuals( mod1, type="deviance" ) 
tempFC$dfbeta<-dfbeta

mod1<-coxph(formula = Surv(tempFC$FCtenure[tempFC$dfbeta<2.5 &tempFC$dfbeta> -2.5], tempFC$FCtenure2[tempFC$dfbeta<2.5 &tempFC$dfbeta> -2.5], tempFC$FCexit[tempFC$dfbeta<2.5 &tempFC$dfbeta> -2.5]) ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempFC[tempFC$dfbeta<2.5 &tempFC$dfbeta> -2.5,],robust=T)

mod2<-coxph(formula = survobjectSRSG_temp ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempSRSG,robust=T)
dfbeta <- residuals( mod2, type="deviance" )
tempSRSG$dfbeta<-dfbeta

mod2<-coxph(formula = Surv(tempSRSG$SRSGtenure[tempSRSG$dfbeta<2.5 &tempSRSG$dfbeta> -2.5], tempSRSG$SRSGtenure2[tempSRSG$dfbeta<2.5 &tempSRSG$dfbeta> -2.5], tempSRSG$SRSGexit[tempSRSG$dfbeta<2.5 &tempSRSG$dfbeta> -2.5]) ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = tempSRSG[tempSRSG$dfbeta<2.5 &tempSRSG$dfbeta> -2.5,],robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]

stargazer(mod1,mod2,digits=2,se=list(robust1,robust2),type="text",df=F,keep.stat=c("n","ll","lr"))


# Table A11

mod1<-glmer(formula = FCexit ~  FCtenure+I(FCtenure^2)+performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+(1|UNacronym), data = dataFC,family="binomial", control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=10000)),nAGQ=0) 
summary(mod1)

mod2<-glmer(formula = SRSGexit ~  SRSGtenure+I(SRSGtenure^2)+performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+(1|UNacronym), data = dataSRSG,family="binomial", control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=10000)),nAGQ=0) 
summary(mod2)

stargazer(mod1,mod2,digits=2,type="text",keep.stat=c("n","ll","aic","bic"))


# Table A12

mod1<-coxph(formula = survobjectFC ~  performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T,method="breslow")

mod2<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataSRSG,robust=T,method="breslow")

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]

stargazer(mod1,mod2,digits=2,se=list(robust1,robust2),type="text",df=F,keep.stat=c("n","ll","lr"))


# Table A13

mod1<-glmer(formula = FCexit ~  FCtenure+I(FCtenure^2)+performance3+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+rotation+(1|UNacronym), data = dataFC,family="binomial",nAGQ=0) 

mod2<-glmer(formula = SRSGexit ~  SRSGtenure+I(SRSGtenure^2)+performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+rotation+missionend+(1|UNacronym), data = dataSRSG,family="binomial",nAGQ=0) 

stargazer(mod1,mod2,digits=2,type="text",keep.stat=c("n","ll","aic","bic"))


# Table A14

mod1<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+log(gle_gdp):log(SRSGtenure)+log(SRSGtenure)+cluster(missionID), data = dataSRSG,robust=T)

mod2<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+SRSGpreviousexperiencedummy:log(SRSGtenure)+log(SRSGtenure)+cluster(missionID), data = dataSRSG,robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]

stargazer(mod1,mod2,digits=2,se=list(robust1,robust2),type="text",df=F,keep.stat=c("n","ll","lr"))


# Table A15

dataFC$performance3_scaled<-scale(dataFC$performance3)
dataSRSG$performance3_scaled<-scale(dataSRSG$performance3)

mod1<-coxph(formula = survobjectFC ~  performance3_scaled+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T)

mod2<-coxph(formula = survobjectFC ~  performance3_scaled+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3_scaled:FC_P5+cluster(missionID), data = dataFC,robust=T)

mod3<-coxph(formula = survobjectFC ~  performance3_scaled+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3_scaled:log(gle_gdp)+cluster(missionID), data = dataFC,robust=T)

mod4<-coxph(formula = survobjectFC ~  performance3_scaled+FC_P5+log(gle_gdp)+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3_scaled:Troops_log+cluster(missionID), data = dataFC,robust=T)

mod5<-coxph(formula = survobjectSRSG ~  performance3_scaled+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataSRSG,robust=T)

mod6<-coxph(formula = survobjectSRSG ~  performance3_scaled+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3_scaled:SRSG_P5+cluster(missionID), data = dataSRSG,robust=T)

mod7<-coxph(formula = survobjectSRSG ~  performance3_scaled+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3_scaled:log(gle_gdp)+cluster(missionID), data = dataSRSG,robust=T)

mod8<-coxph(formula = survobjectSRSG ~  performance3_scaled+SRSG_P5+log(gle_gdp)+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+performance3_scaled:Troops_log+cluster(missionID), data = dataSRSG,robust=T)

robust1     <- summary(mod1)$coefficients[,4]
robust2     <- summary(mod2)$coefficients[,4]
robust3     <- summary(mod3)$coefficients[,4]
robust4     <- summary(mod4)$coefficients[,4]
robust5     <- summary(mod5)$coefficients[,4]
robust6     <- summary(mod6)$coefficients[,4]
robust7     <- summary(mod7)$coefficients[,4]
robust8     <- summary(mod8)$coefficients[,4]

stargazer(mod1,mod2,mod3,mod4,mod5,mod6,mod7,mod8,digits=3,se=list(robust1,robust2,robust3,robust4,robust5,robust6,robust7,robust8),type="text",df=F,keep.stat=c("n","ll","lr"))


# Figures A1 and A2 (maps) not replicated 

# Figure A3, 3 panels

par(mfrow=c(1,2)) 
hist(dataFC$performance3, main="FCs",col="grey", font.main=1, cex.main=.8, xlab="Performance index")
hist(dataSRSG$performance3, main="SRSGs",col="grey", font.main=1, cex.main=.8, xlab="Performance index")
mtext("Performance index", side = 3, line = -1.5, outer = TRUE, font=2, cex=1.2)

par(mfrow=c(1,2)) 
hist(dataFC$Troops, main="FCs",col="grey", font.main=1, cex.main=.8, xlab="Mission troop contribution")
hist(dataSRSG$Troops, main="SRSGs",col="grey", font.main=1, cex.main=.8, xlab="Mission troop contribution")
mtext("Mission troop contribution", side = 3, line = -1.5, outer = TRUE, font=2, cex=1.2)

par(mfrow=c(1,2)) 
hist(dataFC$FC_P5, main="FCs",col="grey", font.main=1, cex.main=.8, xlab="P5 country")
hist(dataSRSG$SRSG_P5, main="SRSGs",col="grey", font.main=1, cex.main=.8, xlab="P5 country")
mtext("P5", side = 3, line = -1.5, outer = TRUE, font=2, cex=1.2)


# Figure A4

mod1<-coxph(formula = survobjectFC ~  performance3+FC_P5+gle_gdp_log+Troops_log+FCpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataFC,robust=T)
mod1
new_data<-data.frame("performance3"=c(62.95,90.35),"Troops_log"=mean(dataFC$Troops_log),FC_P5=0, gle_gdp_log=mean(na.omit(dataFC$gle_gdp_log)),FCpreviousexperiencedummy=0,Chapter.VII=0,Troops_log=na.omit(mean(dataFC$Troops_log)),missionend=0,acting=0)
survfit(mod1,newdata=new_data,conf.int=.95)
par(mfrow=c(1,1))
plot(survfit(mod1,newdata=new_data,conf.int=.90),conf.int=F,lty=c(1,3), xmax=30,lwd=2,main="FCs",xlab="Tenure (months)",ylab="Survival")
legend("topright",lty=c(1,3),legend=c("Low-performing missions","High-performing missions"),lwd=2,cex=.7)

mod5<-coxph(formula = survobjectSRSG ~  performance3+SRSG_P5+gle_gdp+Troops_log+SRSGpreviousexperiencedummy+Chapter.VII+missionend+acting+cluster(missionID), data = dataSRSG,robust=T)
new_data<-data.frame("performance3"=c(49,91.3),"Troops_log"=mean(dataSRSG$Troops_log),SRSG_P5=0, gle_gdp=mean(na.omit(dataSRSG$gle_gdp)),SRSGpreviousexperiencedummy=0,Chapter.VII=0,Troops_log=na.omit(mean(dataSRSG$Troops_log)),missionend=0,acting=0)
survfit(mod5,newdata=new_data,conf.int=.95)
plot(survfit(mod5,newdata=new_data,conf.int=.95),conf.int=F,lty=c(1,3), xmax=30,lwd=2,main="SRSGs",xlab="Tenure (months)",ylab="Survival")
legend("topright",lty=c(1,3),legend=c("Low-performing missions","High-performing missions"),lwd=2,cex=.7)


# Figure A5

# Figure A5: Illustration of performance measures (UNMISS)

par(mfrow=c(6,1))
temp<-dataFC[dataFC$UNacronym=="UNMISS",]
temp$mission_tenure<-temp$mdate-min(temp$mdate)
par(mar=c(1,1,1,1))
plot(1, type="n", xlim=c(min(temp$mission_tenure), max(temp$mission_tenure)), ylim=c(0, 5000), main="Battle-related fatalities",xlab="Mission duration (months)",ylab="Cumulative battle deaths")
for (i in unique(temp$ForceCommander)){
  temp2<-temp[temp$ForceCommander==i,]
  lines(temp2$mission_tenure,temp2$CumulativeFatalities,type="l",lwd=2)
  points(min(temp2$mission_tenure),min(temp2$CumulativeFatalities),pch=16)
}
abline(v=temp$mission_tenure[temp$changeFC==1]-1,col="red",lwd=2)

plot(1, type="n", xlim=c(min(temp$mission_tenure), max(temp$mission_tenure)), ylim=c(0, 10),main="Attacks against mission",xlab="Mission duration (months)",ylab="Cumulative attacks")
for (i in unique(temp$ForceCommander)){
  temp2<-temp[temp$ForceCommander==i,]
  lines(temp2$mission_tenure,temp2$CumulativeMalicious,type="l",lwd=2)
  points(min(temp2$mission_tenure),min(temp2$CumulativeMalicious),pch=16)
}
abline(v=temp$mission_tenure[temp$changeFC==1]-1,col="red",lwd=2)

plot(1, type="n", xlim=c(min(temp$mission_tenure), max(temp$mission_tenure)), ylim=c(0, 1000),main="Fatalities in one-sided violence",xlab="Mission duration (months)",ylab="Cumulative attacks")
for (i in unique(temp$ForceCommander)){
  temp2<-temp[temp$ForceCommander==i,]
  lines(temp2$mission_tenure,temp2$cumulative_onesided_best,type="l",lwd=2)
  points(min(temp2$mission_tenure),min(temp2$cumulative_onesided_best),pch=16)
}

abline(v=temp$mission_tenure[temp$changeFC==1]-1,col="red",lwd=2)
dataFC$cumul
plot(1, type="n",  xlim=c(min(temp$mission_tenure), max(temp$mission_tenure)), ylim=c(0, 100),main="Mission performance index",xlab="Mission duration (months)",ylab="Cumulative events")
for (i in unique(temp$ForceCommander)){
  temp2<-temp[temp$ForceCommander==i,]
  lines(temp2$mission_tenure,temp2$performance3,type="l",lwd=3, col="brown")
  points(min(temp2$mission_tenure),temp2$performance3[1],pch=16)
}
abline(v=temp$mission_tenure[temp$changeFC==1]-1,col="red",lwd=2)


plot(1, type="n",  xlim=c(min(temp$mission_tenure), max(temp$mission_tenure)), ylim=c(-2, 4),main="3-month good performance index",xlab="Mission duration (months)",ylab="Cumulative events")
for (i in unique(temp$ForceCommander)){
  temp2<-temp[temp$ForceCommander==i,]
  lines(temp2$mission_tenure,temp2$cum_good_combined_new_3,type="l",lwd=3, col="brown")
  points(min(temp2$mission_tenure),temp2$cum_good_combined_new_3[1],pch=16)
}
abline(v=temp$mission_tenure[temp$changeFC==1]-1,col="red",lwd=2)

plot(1, type="n",  xlim=c(min(temp$mission_tenure), max(temp$mission_tenure)), ylim=c(-2, 4),main="3-month bad performance index",xlab="Mission duration (months)",ylab="Cumulative events")
for (i in unique(temp$ForceCommander)){
  temp2<-temp[temp$ForceCommander==i,]
  lines(temp2$mission_tenure,temp2$cum_bad_combined_new_3,type="l",lwd=3, col="brown")
  points(min(temp2$mission_tenure),temp2$cum_bad_combined_new_3[1],pch=16)
}
abline(v=temp$mission_tenure[temp$changeFC==1]-1,col="red",lwd=2)


